/*
 * StickNormalAnisoHeatConductionMaterial.C
 *
 *  Created on: 2017年7月5日
 *      Author: liuxiao
 */
/****************************************************************/
/* MOOSE - Multiphysics Object Oriented Simulation Environment  */
/*                                                              */
/*          All contents are licensed under LGPL V2.1           */
/*             See LICENSE for full restrictions                */
/****************************************************************/

#include "StickNormalAnisoHeatConductionMaterial.h"
#include "Function.h"
#include "MooseMesh.h"
#include "Assembly.h"

#include "libmesh/quadrature.h"
registerMooseObject("RoshanApp", StickNormalAnisoHeatConductionMaterial);
template<>
InputParameters validParams<StickNormalAnisoHeatConductionMaterial>()
{
	InputParameters params = validParams<Material>();

	params.addCoupledVar("temperature", "Coupled Temperature");
	params.addCoupledVar("Pe_normal_x", "Coupled Pe_normal_x");
	params.addCoupledVar("Pe_normal_y", "Coupled Pe_normal_y");
	params.addCoupledVar("Pe_normal_z", "Coupled Pe_normal_z");
	params.addCoupledVar("normal_x", "x-component of the local normal");
	params.addCoupledVar("normal_y", "y-component of the local normal");
	params.addCoupledVar("normal_z", "z-component of the local normal");

	params.set<std::vector<VariableName> >("normal_x") = std::vector<VariableName>(1, "nodal_normal_x");
	params.set<std::vector<VariableName> >("normal_y") = std::vector<VariableName>(1, "nodal_normal_y");
	params.set<std::vector<VariableName> >("normal_z") = std::vector<VariableName>(1, "nodal_normal_z");

  params.addParam<Real>("thermal_conductivity_1", "The thermal conductivity in the 1st main direction");
  params.addParam<Real>("thermal_conductivity_2", "The thermal conductivity in the 2nd direction");
  params.addParam<Real>("thermal_conductivity_3", "The thermal conductivity in the 3rd direction");
  params.addParam<Real>("specific_heat", "The specific heat value");
  params.addParam<Real>("density", "The _density value");


  return params;
}



StickNormalAnisoHeatConductionMaterial::StickNormalAnisoHeatConductionMaterial(const InputParameters & parameters) :
  Material(parameters),

  _has_temperature(isCoupled("temperature")),
  _temperature(_has_temperature ? coupledValue("temperature") : _zero),
  _Pe_normal_x(coupledValue("Pe_normal_x") ),
  _Pe_normal_y(coupledValue("Pe_normal_y") ),
  _Pe_normal_z(coupledValue("Pe_normal_z") ),
  _normal_x(coupledValueOld("normal_x")),
  _normal_y(coupledValueOld("normal_y")),
  _normal_z(coupledValueOld("normal_z")),
  _my_thermal_conductivity_1(getParam<Real>("thermal_conductivity_1")),
  _my_thermal_conductivity_2(getParam<Real>("thermal_conductivity_2")),
  _my_thermal_conductivity_3(getParam<Real>("thermal_conductivity_3")),

  _my_specific_heat(isParamValid("specific_heat") ? getParam<Real>("specific_heat") : 0),
  _my_density(isParamValid("density") ? getParam<Real>("density") : 1),



//  _translate_conductivity(declareProperty<RankTwoTensor>("translate_conductivity")),

  _direction_1(declareProperty<Point>("direction_1")),
  _direction_2(declareProperty<Point>("direction_2")),
  _direction_3(declareProperty<Point>("direction_3")),
  _direction_main_x(declareProperty<Real>("direction_main_x")),
  _direction_main_y(declareProperty<Real>("direction_main_y")),
  _direction_main_z(declareProperty<Real>("direction_main_z")),
  _direction_2_x(declareProperty<Real>("direction_2_x")),
  _direction_2_y(declareProperty<Real>("direction_2_y")),
  _direction_2_z(declareProperty<Real>("direction_2_z")),
  _direction_3_x(declareProperty<Real>("direction_3_x")),
  _direction_3_y(declareProperty<Real>("direction_3_y")),
  _direction_3_z(declareProperty<Real>("direction_3_z")),

  _thermal_conductivity_1(declareProperty<Real>("thermal_conductivity_1")),
  _thermal_conductivity_2(declareProperty<Real>("thermal_conductivity_2")),
  _thermal_conductivity_3(declareProperty<Real>("thermal_conductivity_3")),

//  _translate_conductivity(declareProperty<RankTwoTensor>("translate_conductivity")),
//  _conductivity_matrix(declareProperty<RankTwoTensor>("conductivity_matrix")),
//  _translate_conductivity_Trans(declareProperty<RankTwoTensor>("translate_conductivity_Trans")),


  _specific_heat(declareProperty<Real>("specific_heat")),
  _density(declareProperty<Real>("density")),
  _conductivity_matrix_xyz(declareProperty<RankTwoTensor>("thermal_conductivity_xyz"))
{

}

Real
StickNormalAnisoHeatConductionMaterial::compute_cos(Point p1,Point p2)
{
	Real cos_value;
	cos_value = p1*p2/p1.norm()/p2.norm();
	return (cos_value);
}

//RankTwoTensor
//FreeAnisoHeatConductionMaterial::translate_conductivity(Real S11, Real S21, Real S31, Real S12, Real S22, Real S32, Real S13, Real S23, Real S33)
//{
//	RankTwoTensor _translate_conductivity;
//	 _translate_conductivity.RankTwoTensor(S11, S21, S31, S12, S22, S32, S13, S23, S33);
//	return _translate_conductivity;
//}
//
//RankTwoTensor
//FreeAnisoHeatConductionMaterial::conductivity_matrix(Real S11, Real S21, Real S31, Real S12, Real S22, Real S32, Real S13, Real S23, Real S33)
//{
//	RankTwoTensor conductivity_matrix;
//	conductivity_matrix.RankTwoTensor(S11, S21, S31, S12, S22, S32, S13, S23, S33);
//	return conductivity_matrix;
//}




void
StickNormalAnisoHeatConductionMaterial::computeProperties()
{
  for (unsigned int qp(0); qp < _qrule->n_points(); ++qp)
 {
   _thermal_conductivity_1[qp] = _my_thermal_conductivity_1;
   _thermal_conductivity_2[qp] = _my_thermal_conductivity_2;
   _thermal_conductivity_3[qp] = _my_thermal_conductivity_3;
   Point nx(1,0,0);
   Point ny(0,1,0);
   Point nz(0,0,1);
   Point pe(_Pe_normal_x[qp],_Pe_normal_y[qp],_Pe_normal_z[qp]);
   Point pl(_normal_x[qp],_normal_y[qp],_normal_z[qp]);
   Point _my_direction_1(0,0,0);
   if(pe.norm()== 0)
   {
	   _my_direction_1 = pl;
   }
   else
   {
	   _my_direction_1 = pe;
   }

//   std::cout<< _my_direction_1 <<std::endl;
   Point _my_direction_2;
   Point _my_direction_3;
   if(_my_direction_1.cross(nx).norm()>1e-6)
   {
    _my_direction_2 = _my_direction_1.cross(nx);
    _my_direction_3 = -_my_direction_1.cross(_my_direction_2);
   }
   else
   {
	  _my_direction_2 = ny;
	  _my_direction_3 = nz;
   }
   _direction_1[qp] = _my_direction_1.unit();
   _direction_2[qp] = _my_direction_2.unit();
   _direction_3[qp] = _my_direction_3.unit();
   _direction_main_x[qp] = _direction_1[qp](0);
   _direction_main_y[qp] = _direction_1[qp](1);
   _direction_main_z[qp] = _direction_1[qp](2);
   _direction_2_x[qp] = _direction_2[qp](0);
   _direction_2_y[qp] = _direction_2[qp](1);
   _direction_2_z[qp] = _direction_2[qp](2);
   _direction_3_x[qp] = _direction_3[qp](0);
   _direction_3_y[qp] = _direction_3[qp](1);
   _direction_3_z[qp] = _direction_3[qp](2);


   Real l1 = compute_cos(_direction_1[qp],nx);//l1 m1 n1 见文献《正交各向异性材料三维热传导问题的有限元列式》，闫相桥等

   Real l2 = compute_cos(_direction_2[qp],nx);
   Real l3 = compute_cos(_direction_3[qp],nx);
   Real m1 = compute_cos(_direction_1[qp],ny);
   Real m2 = compute_cos(_direction_2[qp],ny);
   Real m3 = compute_cos(_direction_3[qp],ny);
   Real n1 = compute_cos(_direction_1[qp],nz);
   Real n2 = compute_cos(_direction_2[qp],nz);
   Real n3 = compute_cos(_direction_3[qp],nz);
//   RankTwoTensor A(_thermal_conductivity_1[qp], _thermal_conductivity_1[qp], _thermal_conductivity_1[qp]);
   RankTwoTensor translate_conductivity(l1,m1,n1,l2,m2,n2,l3,m3,n3);
   RankTwoTensor conductivity_matrix(_thermal_conductivity_1[qp],0,0,0,_thermal_conductivity_2[qp],0,0,0,_thermal_conductivity_3[qp]);
   RankTwoTensor translate_conductivity_Ttrans(l1,l2,l3,m1,m2,m3,n1,n2,n3);
   RankTwoTensor conductivity_matrix_xyz = translate_conductivity*conductivity_matrix*translate_conductivity_Ttrans;
   _conductivity_matrix_xyz[qp] = conductivity_matrix_xyz;
//   std::cout<<translate_conductivity_Ttrans(0,0)<<","<<translate_conductivity_Ttrans(0,1)<<","<<translate_conductivity_Ttrans(0,2)<<std::endl;
//   std::cout<<translate_conductivity_Ttrans(1,0)<<","<<translate_conductivity_Ttrans(1,1)<<","<<translate_conductivity_Ttrans(1,2)<<std::endl;
//   std::cout<<translate_conductivity_Ttrans(2,0)<<","<<translate_conductivity_Ttrans(2,1)<<","<<translate_conductivity_Ttrans(2,2)<<std::endl;


   _specific_heat[qp] = _my_specific_heat;
   _density[qp] = _my_density;
 }
}





